
import sys
import indexing.featureindex
from collections import deque
import utils.query
from utils.common import printUsage

MAX_READ = 200

def calculateCoverageGenerator(query, index, quitEarlyIfPossible=True):
    """
    Counts the number of mapped bases in the specified summary file
    in the query range.
    """
    (qchr, qmin, qmax) = utils.query.toTuple(query)
    features = indexing.featureindex.getFeatures(index, (qchr, qmin, qmax))
    
    vals = deque([0 for x in xrange(MAX_READ)])
    valsMin = qmin - MAX_READ
    
    cumulativeVal = 0
    for feature in features:
        spos = feature.spos
        epos = feature.epos
        # about to add weight to [spos, epos], so positions left of spos
        # are done.  while valsMin < spos, return a single value and update.
        while valsMin < spos:
            cumulativeVal += vals.popleft()
            if valsMin >= qmin and valsMin <= qmax:
                yield (valsMin, cumulativeVal)
            vals.append(0)
            valsMin += 1
        vals[spos-valsMin] += 1
        vals[epos+1-valsMin] -= 1

    # return remaining non-zero values
    while valsMin <= qmax and len(vals) > 0:
        cumulativeVal += vals.popleft()
        if valsMin>=qmin:
            yield (valsMin, cumulativeVal)
        valsMin += 1
    
    # return zeroes until qmax if that was requested
    if valsMin <= qmax and not quitEarlyIfPossible:
        if valsMin >= qmin:
            yield (valsMin, cumulativeVal)
        valsMin += 1

def getCoverage(query, index, resultArray=None):
    (qchr, qmin, qmax) = utils.query.toTuple(query)
    qlength = qmax - qmin + 1
    if not resultArray:
        resultArray = [0 for x in xrange(qmin, qmax+1)]
    if len(resultArray) < qlength:
        raise Exception("result array not large enough: len=" + len(resultArray) + " vs " + query)
    for x in xrange(qmin, qmax+1):
        resultArray[x-qmin] = 0
    for (x, cov) in calculateCoverageGenerator(query, index):
        if x>=qmin and x<=qmax: resultArray[x-qmin] = cov
    return resultArray


usageString = """
Usage: python %s INDEX_FILE QUERY
Print the number of bases mapped to each position in the query range.
"""

if __name__=="__main__":
    printUsage(sys.argv, usageString % sys.argv[0], 3)
    indexFilename = sys.argv[1]
    query = utils.query.Query(sys.argv[2])
    coverage = calculateCoverageGenerator(query, indexFilename)
    for (x, cov) in coverage:
        print "%s\t%s\t%s" % (query.feature, x, cov)

